roey.angel@bc.cas.cz

Beta diversity analysis

This analysis explores the alpha-diversity ditribution patters in the different samples, based on the DADA2-produced sequences.

Setting general parameters:

set.seed(1000)
min_lib_size <- 5000
data_path <- "./DADA2_pseudo/"
Ps_file <- "TeaTime4Schools_filt3_wTree.Rds"
Proj_name <- "TeaTime4Schools"

Load phyloseq object

This phyloseq object was created in 06_Phylogentic_analysis.html. The Ps_obj_filt object excludes contaminants and all sequences classified as eukaryota, chloroplast, mitochondria or unknown but still includes taxa with low prevalence

readRDS(file = paste0(data_path, Ps_file)) %>%
  subset_samples(., Field != "Unburied") %>% # drop unburied samples
  prune_samples(sample_sums(.) > min_lib_size, .) %>% # remove samples  < min_lib_size
  filter_taxa(., function(x) sum(x) > 0, TRUE) -> # drop taxa with 0 abundance
  Ps_obj_filt
qplot(rowSums(otu_table(Ps_obj_filt)), geom = "histogram") + 
  xlab("Library size")

qplot(log10(rowSums(otu_table(Ps_obj_filt)))) +
  xlab("Logged library size")

Standardize abundances to the median sequencing depth (and convert to proportion)

adonis(
  otu_table(Ps_obj_filt) ~ Lib.size,
  data =
    as(sample_data(Ps_obj_filt), "data.frame"),
  method = "horn",
  permutations = 999
)
## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt) ~ Lib.size, data = as(sample_data(Ps_obj_filt),      "data.frame"), permutations = 999, method = "horn") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Lib.size    1     6.080  6.0804  19.569 0.14122  0.001 ***
## Residuals 119    36.975  0.3107         0.85878           
## Total     120    43.056                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Ps_obj_filt %>%
  otu_table(.) %>%
  as(., "matrix") %>%
  rowSums() %>% 
  median() ->
  total
standf = function(x, t = total) round(t * (x / sum(x)))
Ps_obj_filt_median <- transform_sample_counts(Ps_obj_filt, standf)  # Standardize abundances to median sequencing depth

Ps_obj_filt_median_rel <- transform_sample_counts(Ps_obj_filt_median, function(x) x / sum(x)) # convert to relative abundance (just incase it's explicitly needed)

sample_data(Ps_obj_filt_median)$Lib.size <- sample_sums(Ps_obj_filt_median)

qplot(rowSums(otu_table(Ps_obj_filt_median)), geom = "histogram") + 
  xlab("Library size")

PlotLibDist(Ps_obj_filt_median)

PlotReadHist(as(otu_table(Ps_obj_filt_median), "matrix"))

notAllZero <- (rowSums(t(otu_table(Ps_obj_filt_median))) > 0)
meanSdPlot(as(log2(t(otu_table(Ps_obj_filt_median))[notAllZero, ] + 1), "matrix"))

Variance partitioning models and ordinations

Partitioning the data using discrete distance

adonis(
  otu_table(Ps_obj_filt_median) ~ Lib.size,
  data =
    as(sample_data(Ps_obj_filt_median), "data.frame"),
  method = "horn",
  permutations = 999
)
## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median) ~ Lib.size, data = as(sample_data(Ps_obj_filt_median),      "data.frame"), permutations = 999, method = "horn") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Lib.size    1     0.285 0.28485 0.79255 0.00662   0.57
## Residuals 119    42.770 0.35941         0.99338       
## Total     120    43.055                 1.00000
adonis(
  otu_table(Ps_obj_filt_median_rel) ~ Field + Sample.type * Season,
  data =
    as(sample_data(Ps_obj_filt_median_rel), "data.frame"),
  method = "horn",
  permutations = 999
)
## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median_rel) ~ Field +      Sample.type * Season, data = as(sample_data(Ps_obj_filt_median_rel),      "data.frame"), permutations = 999, method = "horn") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                     Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field                2     1.158  0.5789   5.305 0.02689  0.001 ***
## Sample.type          2    16.495  8.2477  75.584 0.38313  0.001 ***
## Season               3     8.674  2.8913  26.497 0.20146  0.001 ***
## Sample.type:Season   6     5.052  0.8419   7.716 0.11733  0.001 ***
## Residuals          107    11.676  0.1091         0.27119           
## Total              120    43.055                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# sample_data(Ps_obj_filt_median)$Type.season <- paste(sample_data(Ps_obj_filt)$Sample.type, sample_data(Ps_obj_filt)$Season) # too many comparisons

# Ps_obj_filt_median_s <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median), TRUE)[1:100]), Ps_obj_filt_median) # for testing

# Compare sample-type pairs
mod_pairwise1 <- PairwiseAdonis(
  otu_table(Ps_obj_filt_median),
  sample_data(Ps_obj_filt_median)$Sample.type,
  sim.function = "vegdist",
  sim.method = "horn",
  p.adjust.m = "BH"
)
print(mod_pairwise1)
##                  pairs total.DF  F.Model        R2 p.value p.adjusted sig
## 1 Green tea vs Rooibos       91 13.90875 0.1338555   0.001      0.001  **
## 2    Green tea vs Soil       73 48.04850 0.4002424   0.001      0.001  **
## 3      Rooibos vs Soil       75 67.00992 0.4752142   0.001      0.001  **
(sig_pairs1 <- as.character(mod_pairwise1$pairs[mod_pairwise1$p.adjusted < 0.05]))
## [1] "Green tea vs Rooibos" "Green tea vs Soil"    "Rooibos vs Soil"
# (meaningful_sig_pairs1 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))

# compare seasons
mod_pairwise2 <- PairwiseAdonis(
  otu_table(Ps_obj_filt_median),
  sample_data(Ps_obj_filt_median)$Season,
  sim.function = "vegdist",
  sim.method = "horn",
  p.adjust.m = "BH"
)
print(mod_pairwise2)
##              pairs total.DF   F.Model         R2 p.value p.adjusted sig
## 1 Winter vs Spring       58  7.531288 0.11670754   0.001     0.0012   *
## 2 Winter vs Summer       59 21.038047 0.26617620   0.001     0.0012   *
## 3 Winter vs Autumn       61 16.210974 0.21271181   0.001     0.0012   *
## 4 Spring vs Summer       58  7.851482 0.12106866   0.001     0.0012   *
## 5 Spring vs Autumn       60  6.103328 0.09374832   0.001     0.0012   *
## 6 Summer vs Autumn       61  5.653014 0.08610441   0.002     0.0020   *
(sig_pairs2 <- as.character(mod_pairwise2$pairs[mod_pairwise2$p.adjusted < 0.05]))
## [1] "Winter vs Spring" "Winter vs Summer" "Winter vs Autumn" "Spring vs Summer"
## [5] "Spring vs Autumn" "Summer vs Autumn"
# (meaningful_sig_pairs2 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))
Calculate ordinations
Ps_obj_ord1 <- ordinate(Ps_obj_filt_median, "CAP", "horn", formula = Ps_obj_filt_median ~  Field + Sample.type * Season)
evals <- eigenvals(Ps_obj_ord1) # /sum(eigenvals(Ps_obj_ord)) * 100

Ps_obj_filt_median %>% 
  plot_ordination(., Ps_obj_ord1, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>% 
  mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>% 
  mutate_at(., "Sample.type", ~fct_relevel(., "Soil", "Green tea", "Rooibos")) %>% 
  dplyr::rename(., `Sample type` = Sample.type) ->
  ord_df

# p_ord.file <- "PCoA_bray"
# svglite(paste0(p_ord.file, ".svg"),
#         width = 10, height = 7.2)

p_ord <- ggplot(ord_df,
             aes(
               x = CAP1,
               y = CAP2,
               shape = Field,
               color = `Sample type`
             )) +
  geom_point(size = 4, alpha = 2 / 3) +
  theme_bw(base_size = 14) +
  scale_colour_manual(values = pal_locuszoom("default")(3)[c(2, 3, 1)]) +
  stat_ellipse(
    aes(x = CAP1, y = CAP2, group = `Sample type`),
    alpha = 0.5,
    type = "norm",
    level = 0.95,
    linetype = 2,
    inherit.aes = FALSE
  ) +
  labs(x = sprintf("CAP1 [%s%%]", round(evals[1], 1)), 
       y = sprintf("CAP2 [%s%%]", round(evals[2], 1))) +
  coord_fixed(sqrt(evals[2] / evals[1])) +
  facet_wrap(~ Season)

# # Now add the environmental variables as arrows
# arrowmat <- scores(Ps_obj_ord1, display = "bp") # bipplot arrows
# # Add labels, make a data.frame
# arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)
# 
# arrowdf %<>%
#   mutate(labels = fct_recode(labels,
#     "Braunerde" = "FieldBraunerde",
#     "Kolluvisol" = "FieldKolluvisol",
#     "Rooibos" = "Sample.typeRooibos",
#     "Soil" = "Sample.typeSoil"
#   ))
# arrowdf %<>% dplyr::slice(., c(1:4))
# 
# # Define the arrow aesthetic mapping
# arrow_map = aes(xend = CAP1, yend = CAP2, x = 0, y = 0, shape = NULL, color = NULL, 
#     label = labels)
# label_map = aes(x = 1.2 * CAP1, y = 1.2 * CAP2, shape = NULL, color = NULL, 
#     label = labels)
# # Make a new graphic
# arrowhead = arrow(length = unit(0.05, "npc"))
# p_ord <- p_ord +
#   geom_segment(
#     arrow_map,
#     size = 0.5,
#     data = arrowdf,
#     color = "gray",
#     arrow = arrowhead
#   ) + 
#   geom_text(label_map, size = 2, data = arrowdf)
print(p_ord)

# dev.off()

# ggsave(
#   paste0(p_ord.file, ".png"),
#   p_ord,
#   device = "png",
#   width = 10,
#   height = 6
# )
# gz(paste0(p_ord.file, ".svg"), paste0(p_ord.file, ".svgz"))

Partitioning the data using phylometric distance

Unifrac_mat <- UniFrac(Ps_obj_filt_median, 
                       weighted = TRUE, 
                       normalized = TRUE, 
                       parallel = TRUE, 
                       fast = TRUE)

adonis(
  Unifrac_mat ~ Lib.size,
  data =
    as(sample_data(Ps_obj_filt_median), "data.frame"),
  permutations = 999
)
## 
## Call:
## adonis(formula = Unifrac_mat ~ Lib.size, data = as(sample_data(Ps_obj_filt_median),      "data.frame"), permutations = 999) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## Lib.size    1    0.0589 0.058917  1.2597 0.01047  0.254
## Residuals 119    5.5658 0.046771         0.98953       
## Total     120    5.6247                  1.00000
adonis(
  Unifrac_mat ~ Field + Sample.type * Season,
  data =
    as(sample_data(Ps_obj_filt_median), "data.frame"),
  permutations = 999
)
## 
## Call:
## adonis(formula = Unifrac_mat ~ Field + Sample.type * Season,      data = as(sample_data(Ps_obj_filt_median), "data.frame"),      permutations = 999) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                     Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field                2    0.0703 0.03515   3.295 0.01250  0.007 ** 
## Sample.type          2    2.6715 1.33576 125.228 0.47496  0.001 ***
## Season               3    1.1516 0.38388  35.989 0.20475  0.001 ***
## Sample.type:Season   6    0.5899 0.09832   9.218 0.10488  0.001 ***
## Residuals          107    1.1413 0.01067         0.20291           
## Total              120    5.6247                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Calculate ordinations
Ps_obj_ord2 <- ordinate(Ps_obj_filt_median, "PCoA", "unifrac", weighted = TRUE, formula = Ps_obj_filt_median ~ Field + Sample.type * Season)
evals <- Ps_obj_ord2$values$Eigenvalues[1:2] # /sum(eigenvals(Ps_obj_ord)) * 100

Ps_obj_filt_median %>% 
  plot_ordination(., Ps_obj_ord2, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>% 
  mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>% 
  mutate_at(., "Sample.type", ~fct_relevel(., "Soil", "Green tea", "Rooibos")) %>% 
  dplyr::rename(., `Sample type` = Sample.type) ->
  ord_df2

# p_ord2.file <- "PCoA_bray"
# svglite(paste0(p_ord2.file, ".svg"),
#         width = 10, height = 7.2)

p_ord2 <- ggplot(ord_df2,
             aes(
               x = Axis.1,
               y = Axis.2,
               shape = Field,
               color = `Sample type`
             )) +
  geom_point(size = 4, alpha = 2 / 3) +
  theme_bw(base_size = 14) +
  scale_colour_manual(values = pal_locuszoom("default")(3)[c(2, 3, 1)]) +
  stat_ellipse(
    aes(x = Axis.1, y = Axis.2, group = `Sample type`),
    alpha = 0.5,
    type = "norm",
    level = 0.95,
    linetype = 2,
    inherit.aes = FALSE
  ) +
  labs(x = sprintf("Axis1 [%s%%]", round(evals[1], 1)), 
       y = sprintf("Axis2 [%s%%]", round(evals[2], 1))) +
  coord_fixed(sqrt(evals[2] / evals[1])) +
  facet_wrap( ~ Season)

# plot_grid(p_ord,p_ord2)
print(p_ord2)

# dev.off()

# ggsave(
#   paste0(p_ord2.file, ".png"),
#   p_ord2,
#   device = "png",
#   width = 10,
#   height = 6
# )
# gz(paste0(p_ord2.file, ".svg"), paste0(p_ord2.file, ".svgz"))

Variance partitioning models and ordinations - soils excluded

Partitioning the data using discrete distance

Ps_obj_filt_median_tea <- subset_samples(Ps_obj_filt_median, Sample.type != "Soil")

adonis(
  otu_table(Ps_obj_filt_median_tea) ~ Lib.size,
  data =
    as(sample_data(Ps_obj_filt_median_tea), "data.frame"),
  method = "horn",
  permutations = 999
)
## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median_tea) ~ Lib.size,      data = as(sample_data(Ps_obj_filt_median_tea), "data.frame"),      permutations = 999, method = "horn") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## Lib.size   1    0.3309 0.33093  1.0424 0.01145  0.416
## Residuals 90   28.5712 0.31746         0.98855       
## Total     91   28.9021                 1.00000
adonis(
  otu_table(Ps_obj_filt_median_tea) ~ Field + Sample.type * Season,
  data =
    as(sample_data(Ps_obj_filt_median_tea), "data.frame"),
  method = "horn",
  permutations = 999
)
## 
## Call:
## adonis(formula = otu_table(Ps_obj_filt_median_tea) ~ Field +      Sample.type * Season, data = as(sample_data(Ps_obj_filt_median_tea),      "data.frame"), permutations = 999, method = "horn") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                    Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field               2    1.3803  0.6901  5.4738 0.04776  0.001 ***
## Sample.type         1    3.8840  3.8840 30.8053 0.13439  0.001 ***
## Season              3   10.4768  3.4923 27.6981 0.36249  0.001 ***
## Sample.type:Season  3    2.8222  0.9407  7.4612 0.09765  0.001 ***
## Residuals          82   10.3388  0.1261         0.35772           
## Total              91   28.9021                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# sample_data(Ps_obj_filt_median)$Type.season <- paste(sample_data(Ps_obj_filt)$Sample.type, sample_data(Ps_obj_filt)$Season) # too many comparisons

# Ps_obj_filt_median_s <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median), TRUE)[1:100]), Ps_obj_filt_median) # for testing

# Compare sample-type pairs
mod_pairwise3 <- PairwiseAdonis(
  otu_table(Ps_obj_filt_median_tea),
  sample_data(Ps_obj_filt_median_tea)$Sample.type,
  sim.function = "vegdist",
  sim.method = "horn",
  p.adjust.m = "BH"
)
print(mod_pairwise3)
##                  pairs total.DF  F.Model        R2 p.value p.adjusted sig
## 1 Green tea vs Rooibos       91 13.90875 0.1338555   0.001      0.001  **
(sig_pairs3 <- as.character(mod_pairwise3$pairs[mod_pairwise3$p.adjusted < 0.05]))
## [1] "Green tea vs Rooibos"
# (meaningful_sig_pairs3 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))

# compare seasons
mod_pairwise4 <- PairwiseAdonis(
  otu_table(Ps_obj_filt_median_tea),
  sample_data(Ps_obj_filt_median_tea)$Season,
  sim.function = "vegdist",
  sim.method = "horn",
  p.adjust.m = "BH"
)
print(mod_pairwise4)
##              pairs total.DF  F.Model        R2 p.value p.adjusted sig
## 1 Winter vs Spring       46 13.19236 0.2267026   0.001      0.001  **
## 2 Winter vs Summer       44 34.98156 0.4485876   0.001      0.001  **
## 3 Winter vs Autumn       47 28.36909 0.3814635   0.001      0.001  **
## 4 Spring vs Summer       43 11.80495 0.2194027   0.001      0.001  **
## 5 Spring vs Autumn       46 10.05915 0.1826972   0.001      0.001  **
## 6 Summer vs Autumn       44  8.91082 0.1716563   0.001      0.001  **
(sig_pairs4 <- as.character(mod_pairwise4$pairs[mod_pairwise4$p.adjusted < 0.05]))
## [1] "Winter vs Spring" "Winter vs Summer" "Winter vs Autumn" "Spring vs Summer"
## [5] "Spring vs Autumn" "Summer vs Autumn"
# (meaningful_sig_pairs2 <- c("Green tea vs Rooibos", "Green tea vs Soil", "Rooibos vs Soil"))
Calculate ordinations
Ps_obj_ord3 <- ordinate(Ps_obj_filt_median_tea, "CAP", "horn", formula = Ps_obj_filt_median_tea ~  Field + Sample.type * Season)
evals <- eigenvals(Ps_obj_ord3) # /sum(eigenvals(Ps_obj_ord)) * 100

Ps_obj_filt_median_tea %>% 
  plot_ordination(., Ps_obj_ord3, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>% 
  mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>% 
  mutate_at(., "Sample.type", ~fct_relevel(., "Green tea", "Rooibos")) %>% 
  dplyr::rename(., `Sample type` = Sample.type) ->
  ord_df3

# p_ord.file <- "PCoA_bray"
# svglite(paste0(p_ord.file, ".svg"),
#         width = 10, height = 7.2)

p_ord3 <- ggplot(ord_df3,
             aes(
               x = CAP1,
               y = CAP2,
               shape = Field,
               color = `Sample type`
             )) +
  geom_point(size = 4, alpha = 2 / 3) +
  theme_bw(base_size = 14) +
  scale_colour_manual(values = pal_locuszoom("default")(3)[c(3, 1)]) +
  stat_ellipse(
    aes(x = CAP1, y = CAP2, group = `Sample type`),
    alpha = 0.5,
    type = "norm",
    level = 0.95,
    linetype = 2,
    inherit.aes = FALSE
  ) +
  labs(x = sprintf("Axis1 [%s%%]", round(evals[1], 1)), 
       y = sprintf("Axis2 [%s%%]", round(evals[2], 1))) +
  coord_fixed(sqrt(evals[2] / evals[1])) +
  facet_wrap(~ Season)

print(p_ord3)

# dev.off()

# ggsave(
#   paste0(p_ord.file, ".png"),
#   p_ord,
#   device = "png",
#   width = 10,
#   height = 6
# )
# gz(paste0(p_ord.file, ".svg"), paste0(p_ord.file, ".svgz"))

Partitioning the data using phylometric distance

Unifrac_mat <- UniFrac(Ps_obj_filt_median_tea, 
                       weighted = TRUE, 
                       normalized = TRUE, 
                       parallel = TRUE, 
                       fast = TRUE)

adonis(
  Unifrac_mat ~ Lib.size,
  data =
    as(sample_data(Ps_obj_filt_median_tea), "data.frame"),
  permutations = 999
)
## 
## Call:
## adonis(formula = Unifrac_mat ~ Lib.size, data = as(sample_data(Ps_obj_filt_median_tea),      "data.frame"), permutations = 999) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)  
## Lib.size   1    0.0431 0.043104  2.0314 0.02207  0.076 .
## Residuals 90    1.9097 0.021219         0.97793         
## Total     91    1.9528                  1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(
  Unifrac_mat ~ Field + Sample.type * Season,
  data =
    as(sample_data(Ps_obj_filt_median_tea), "data.frame"),
  permutations = 999
)
## 
## Call:
## adonis(formula = Unifrac_mat ~ Field + Sample.type * Season,      data = as(sample_data(Ps_obj_filt_median_tea), "data.frame"),      permutations = 999) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                    Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)    
## Field               2   0.04872 0.024358   3.170 0.02495  0.008 ** 
## Sample.type         1   0.16444 0.164441  21.403 0.08421  0.001 ***
## Season              3   0.94795 0.315983  41.127 0.48543  0.001 ***
## Sample.type:Season  3   0.16168 0.053893   7.015 0.08279  0.001 ***
## Residuals          82   0.63001 0.007683         0.32262           
## Total              91   1.95279                  1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Calculate ordinations
Ps_obj_ord4 <- ordinate(Ps_obj_filt_median_tea, "PCoA", "unifrac", weighted = TRUE, formula = Ps_obj_filt_median_tea ~ Field + Sample.type * Season)
evals <- Ps_obj_ord4$values$Eigenvalues[1:2] # /sum(eigenvals(Ps_obj_ord)) * 100

Ps_obj_filt_median_tea %>% 
  plot_ordination(., Ps_obj_ord4, type = "samples", shape = "Field", color = "Sample.type", justDF = TRUE) %>% 
  mutate_at(., "Season", ~fct_relevel(., "Winter", "Spring", "Summer", "Autumn")) %>% 
  mutate_at(., "Sample.type", ~fct_relevel(., "Green tea", "Rooibos")) %>% 
  dplyr::rename(., `Sample type` = Sample.type) ->
  ord_df4

# p_ord.file <- "PCoA_bray"
# svglite(paste0(p_ord.file, ".svg"),
#         width = 10, height = 7.2)

p_ord4 <- ggplot(ord_df4,
             aes(
               x = Axis.1,
               y = Axis.2,
               shape = Field,
               color = `Sample type`
             )) +
  geom_point(size = 4, alpha = 2 / 3) +
  theme_bw(base_size = 14) +
  scale_colour_manual(values = pal_locuszoom("default")(3)[c(3, 1)]) +
  stat_ellipse(
    aes(x = Axis.1, y = Axis.2, group = `Sample type`),
    alpha = 0.5,
    type = "norm",
    level = 0.95,
    linetype = 2,
    inherit.aes = FALSE
  ) +
  labs(x = sprintf("Axis1 [%s%%]", round(evals[1], 1)), 
       y = sprintf("Axis2 [%s%%]", round(evals[2], 1))) +
  coord_fixed(sqrt(evals[2] / evals[1])) +
  facet_wrap( ~ Season)

# plot_grid(p_ord,p_ord)
print(p_ord4)

# dev.off()

# ggsave(
#   paste0(p_ord.file, ".png"),
#   p_ord,
#   device = "png",
#   width = 10,
#   height = 6
# )
# gz(paste0(p_ord.file, ".svg"), paste0(p_ord.file, ".svgz"))

Test differences between samples on the phylum level

Taxa_tests_phylum1 <- STAMPR(Ps_obj_filt_median, vars2test = "Sample.type", sig_pairs = sig_pairs1, outputfile = paste0(Proj_name, "_Sample.type"))
plotSTAMPR(Taxa_tests_phylum1, pair = "Green tea vs Rooibos")

plotSTAMPR(Taxa_tests_phylum1, pair = "Green tea vs Soil")

plotSTAMPR(Taxa_tests_phylum1, pair = "Rooibos vs Soil")

Taxa_tests_order1 <- STAMPR(Ps_obj_filt_median, vars2test = "Sample.type", rank = "Order", sig_pairs = sig_pairs1, outputfile = paste0(Proj_name, "_Sample.type"))
plotSTAMPR(Taxa_tests_order1, pair = "Green tea vs Rooibos", tax_level = "Order")

plotSTAMPR(Taxa_tests_order1, pair = "Green tea vs Soil", tax_level = "Order")

plotSTAMPR(Taxa_tests_order1, pair = "Rooibos vs Soil", tax_level = "Order")

Taxa_tests_phylum2 <- STAMPR(Ps_obj_filt_median, vars2test = "Season", sig_pairs = sig_pairs2, outputfile = paste0(Proj_name, "_Season"))
plotSTAMPR(Taxa_tests_phylum2, pair = "Winter vs Summer")

Taxa_tests_order2 <- STAMPR(Ps_obj_filt_median, vars2test = "Season", rank = "Order", sig_pairs = sig_pairs2, outputfile = paste0(Proj_name, "_Season"))
plotSTAMPR(Taxa_tests_order2, pair = "Winter vs Summer", tax_level = "Order")

Taxa_tests_phylum4 <- STAMPR(Ps_obj_filt_median_tea, vars2test = "Season", sig_pairs = sig_pairs4, outputfile = paste0(Proj_name, "Teabags_Season"))
plotSTAMPR(Taxa_tests_phylum4, pair = "Winter vs Summer")

Taxa_tests_order4 <- STAMPR(Ps_obj_filt_median_tea, vars2test = "Season", rank = "Order", sig_pairs = sig_pairs4, outputfile = paste0(Proj_name, "Teabags_Season"))
plotSTAMPR(Taxa_tests_order4, pair = "Winter vs Summer", tax_level = "Order")

Differential abundance models

Tag rare phyla (for plotting purposes only)

Ps_obj_filt_media_glom <- tax_glom(Ps_obj_filt_median, 
                             "Phylum", 
                             NArm = TRUE) # glomerate to the phylum level
Ps_obj_filt_media_glom_rel <- transform_sample_counts(Ps_obj_filt_media_glom, function(x) x / sum(x)) # transform to rel. ab.
Ps_obj_filt_media_glom_rel_DF <- psmelt(Ps_obj_filt_media_glom_rel) # generate a df
Ps_obj_filt_media_glom_rel_DF$Phylum %<>% as.character() # factor to char

# group dataframe by Phylum, calculate median rel. abundance
Ps_obj_filt_media_glom_rel_DF %>%
  group_by(Phylum) %>%
  summarise(median = median(Abundance)) ->
  medians

# find Phyla whose median rel. abund. is less than 0.5%
Rare_phyla <- medians[medians$median <= 0.005, ]$Phylum

# change their name to "Rare"
Ps_obj_filt_media_glom_rel_DF[Ps_obj_filt_media_glom_rel_DF$Phylum %in% Rare_phyla, ]$Phylum <- 'Rare'

# re-group
Ps_obj_filt_media_glom_rel_DF %>%
  group_by(Phylum) %>%
  summarise(Abundance = sum(Abundance)) %>% 
  arrange(desc(Abundance)) -> Taxa_rank

Detect differentially abundant OTUs using ALDEx2 (Fernandes et al. 2013)

Sample type differences

significance = 0.05

# run full model 
data2test <- t(otu_table(Ps_obj_filt_median))
# comparison <- as.character(get_variable(Ps_obj_filt_median, "Sample.type"))
#ALDEx_full <- aldex.clr(data2test, comparison, mc.samples = 128, denom = "iqlr", verbose = TRUE, useMC = TRUE) # iqlr for slight assymetry in composition
#ALDEx_full_glm <- aldex.glm(ALDEx_full, comparison, useMC = TRUE) # for more than two conditions
#sig_taxa <- rownames(ALDEx_full_glm)[ALDEx_full_glm$glm.eBH < 0.05] # save names of taxa that are significant under the full model
#write.csv(sig_taxa, file = "Aldex_full_significant_taxa.csv")

# Pairwise comparisons
# 
ALDEx_comparisons <- list()
ALDEx_comparisons$Comparisons <- sig_pairs1

# Ps_obj_filt_median_subset <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median), TRUE)[1:100]), Ps_obj_filt_median)

for (j in seq(1, length(ALDEx_comparisons$Comparisons))) {
  # print(j)
  ALDEx_comparisons$Comparisons[j] %>% 
    str_split(., " vs ", simplify = FALSE) %>% 
    unlist() ->
    comparison_string
  Ps_obj_filt_median %>%
    subset_samples(Sample.type == comparison_string[1] | Sample.type == comparison_string[2]) ->
    # subset_samples(Treatment == comparison_string[2] | Treatment == comparison_string[4]) %>%
    # subset_samples(Year == "2016" | Year == "2017" | Year == "2018") ->
    Ps_obj_filt_median_pairwise

#  Remove species with prevalence < 10%
  Ps_obj_filt_median_pairwise_s <- DropRareSpecies(Ps_obj = Ps_obj_filt_median_pairwise, prevalence = 0.1)
    
  sample_data(Ps_obj_filt_median_pairwise_s)$Sample.type %<>% fct_relevel(., str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1])
  
  # make Joint.sample.name for matching OTUs between compared samples (for GGPlotTopOTUs)
  suppressWarnings(
  sample_data(Ps_obj_filt_median_pairwise_s) %<>%
    as(., "data.frame") %>% 
    rownames_to_column() %>% 
    mutate_if(is.factor, as.character) %>% 
    mutate(Joint.sample.name  = paste0(.$Sample.type, "_", .$Field, "_", .$Season, "_", .$Replicate)) %>% 
    mutate_at("Joint.sample.name", ~str_replace_all(., paste0(unique(comparison_string), collapse = "|"), "")) %>% # remove the levels participating in the comparison from the names
    mutate_at("Joint.sample.name", ~str_replace_all(., "^_", "")) %>% # remove first "_"
    column_to_rownames()
  )
    
  ALDEx2plot_pairwise <- CalcALDEx(
    physeq_obj = Ps_obj_filt_median_pairwise_s,
    vars2test = "Sample.type",
    rare_phyla = Rare_phyla,
    sig_level = significance,
    LFC = 0
    )
  ALDEx_comparisons$Results[[j]] <- ALDEx2plot_pairwise # store results
  ALDEx_comparisons$Results[[j]]$Var1 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1]
  ALDEx_comparisons$Results[[j]]$Var2 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[2]
  ALDEX_summary <- tibble(Label = c(paste0("⬆", sum( ALDEx_comparisons$Results[[j]]$effect > 0 &  ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " ⬇", sum( ALDEx_comparisons$Results[[j]]$effect < 0 &  ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " (", nrow( ALDEx_comparisons$Results[[j]]), ")")))
  
  ALDEx2plot_pairwise %>%
    filter(Significance == "Pass") %>%
    dplyr::select(OTU, baseMean, effect, Phylum, Class, Order, Family, Genus) %>%
    arrange(desc(abs(effect))) ->
    ALDEx2plot_pairwise_results
    
  # print(ALDEx2plot_pairwise_results %>% 
  #   kable(., digits = c(2), caption = "Significantly different taxa:") %>%
  #   kable_styling(
  #     bootstrap_options = c("striped", "hover", "condensed", "responsive"),
  #     full_width = F
  #   ))
  
  write.csv(ALDEx2plot_pairwise_results, file = paste0("Aldex", "_", paste0(comparison_string, collapse = "_"), ".csv"))
  
  # Plot ALDEX plot
  p1 <- GGPlotALDExTax(ALDEx2plot_pairwise, OTU_labels = FALSE, sig_level = significance) +
    # ggtitle(ALDEx_comparisons$Comparisons[j]) +
    geom_text(
    data    = ALDEX_summary,
    mapping = aes(x = Inf, y = Inf, label = Label),
    hjust   = 1.1,
    vjust   = 1.6
  ) 
  p1 <- p1 + labs(title = ALDEx_comparisons$Comparisons[j])
  print(p1)
  
  # Plot OTU plots
  # GGPlotOTU(Ps_obj_filt_subset_pairwise_s, vars2test = "Spill.Treatment", "Seq_8")
  p2 <- GGPlotTopOTUs(
    physeq_obj = Ps_obj_filt_median_pairwise_s,
    vars2test = "Sample.type",
    ALDEx_obj = ALDEx2plot_pairwise,
    rank_by = "effect",
    Ntop = 12
  )
  print(p2)
}

phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 92 samples ] sample_data() Sample Data: [ 92 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3042 taxa and 92 samples ] sample_data() Sample Data: [ 92 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3042 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3042 tips and 3040 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 74 samples ] sample_data() Sample Data: [ 74 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3217 taxa and 74 samples ] sample_data() Sample Data: [ 74 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3217 tips and 3215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 76 samples ] sample_data() Sample Data: [ 76 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3411 taxa and 76 samples ] sample_data() Sample Data: [ 76 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3411 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3411 tips and 3409 internal nodes ]

filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU[common_taxa <- filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU %in% filter(ALDEx_comparisons$Results[[2]], Significance == "Pass")$OTU]

[1] "Seq_64" "Seq_475" "Seq_631" "Seq_620" "Seq_120" "Seq_112" "Seq_797" [8] "Seq_169" "Seq_37" "Seq_124" "Seq_170" "Seq_1032" "Seq_1382" "Seq_403" [15] "Seq_1094" "Seq_763" "Seq_644" "Seq_855" "Seq_947" "Seq_922" "Seq_507" [22] "Seq_75" "Seq_1003" "Seq_659" "Seq_116" "Seq_698" "Seq_323" "Seq_553" [29] "Seq_313" "Seq_733" "Seq_152" "Seq_968" "Seq_1165" "Seq_838" "Seq_25"
[36] "Seq_2055" "Seq_964" "Seq_434" "Seq_469" "Seq_436" "Seq_65" "Seq_83"
[43] "Seq_1042" "Seq_800" "Seq_488" "Seq_237" "Seq_446" "Seq_264" "Seq_535" [50] "Seq_101" "Seq_665" "Seq_268" "Seq_192" "Seq_940" "Seq_281" "Seq_1305" [57] "Seq_166" "Seq_364" "Seq_2092" "Seq_847" "Seq_868" "Seq_175" "Seq_76"
[64] "Seq_68" "Seq_316" "Seq_122" "Seq_496" "Seq_23" "Seq_195" "Seq_17"
[71] "Seq_15" "Seq_202" "Seq_348" "Seq_270" "Seq_431" "Seq_154" "Seq_61"
[78] "Seq_1232" "Seq_709" "Seq_1048" "Seq_214" "Seq_358" "Seq_1128" "Seq_1001" [85] "Seq_854" "Seq_159" "Seq_42" "Seq_439" "Seq_419" "Seq_33" "Seq_380" [92] "Seq_1058" "Seq_1172" "Seq_367" "Seq_397" "Seq_126" "Seq_132" "Seq_483" [99] "Seq_52" "Seq_306" "Seq_115" "Seq_420" "Seq_31" "Seq_41" "Seq_251" [106] "Seq_12" "Seq_5" "Seq_56" "Seq_2" "Seq_109" "Seq_534" "Seq_131" [113] "Seq_664" "Seq_452" "Seq_208" "Seq_227" "Seq_1404" "Seq_786" "Seq_287" [120] "Seq_235" "Seq_646" "Seq_1057" "Seq_801" "Seq_223" "Seq_635" "Seq_310" [127] "Seq_489" "Seq_318" "Seq_48" "Seq_501" "Seq_438" "Seq_95" "Seq_1738" [134] "Seq_1756" "Seq_2179" "Seq_645" "Seq_1294" "Seq_984" "Seq_788" "Seq_1470" [141] "Seq_228" "Seq_111" "Seq_372" "Seq_253" "Seq_10" "Seq_58" "Seq_1685" [148] "Seq_1186" "Seq_447" "Seq_1099" "Seq_130" "Seq_11" "Seq_578" "Seq_458" [155] "Seq_595" "Seq_565" "Seq_6" "Seq_840" "Seq_2786"

Seasonal differences Comparing the seasonsal differences in teabag samples only

significance = 0.05

data2test <- t(otu_table(Ps_obj_filt_median_tea))

# Pairwise comparisons
ALDEx_comparisons <- list()
ALDEx_comparisons$Comparisons <- sig_pairs2

# Ps_obj_filt_median_tea_subset <- prune_taxa(names(sort(taxa_sums(Ps_obj_filt_median_tea), TRUE)[1:100]), Ps_obj_filt_median_tea)

for (j in seq(1, length(ALDEx_comparisons$Comparisons))) {
  # print(j)
  ALDEx_comparisons$Comparisons[j] %>% 
    str_split(., " vs ", simplify = FALSE) %>% 
    unlist() ->
    comparison_string
  Ps_obj_filt_median_tea %>%
    subset_samples(Season == comparison_string[1] | Season == comparison_string[2]) ->
    # subset_samples(Treatment == comparison_string[2] | Treatment == comparison_string[4]) %>%
    # subset_samples(Year == "2016" | Year == "2017" | Year == "2018") ->
    Ps_obj_filt_median_tea_pairwise

#  Remove species with prevalence < 10%
  Ps_obj_filt_median_tea_pairwise_s <- DropRareSpecies(Ps_obj = Ps_obj_filt_median_tea_pairwise, prevalence = 0.1)
    
  sample_data(Ps_obj_filt_median_tea_pairwise_s)$Season %<>% fct_relevel(., str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1])
  
  # make Joint.sample.name for matching OTUs between compared samples (for GGPlotTopOTUs)
  suppressWarnings(
  sample_data(Ps_obj_filt_median_tea_pairwise_s) %<>% 
    as(., "data.frame") %>% 
    rownames_to_column() %>% 
    mutate_if(is.factor, as.character) %>% 
    mutate(Joint.sample.name  = paste0(.$Sample.type, "_", .$Field, "_", .$Season, "_", .$Replicate)) %>% 
    mutate_at("Joint.sample.name", ~str_replace_all(., paste0(unique(comparison_string), collapse = "|"), "")) %>% # remove the levels participating in the comparison from the names
    mutate_at("Joint.sample.name", ~str_replace_all(., "^_", "")) %>% # remove first "_"
    mutate_at("Joint.sample.name", ~str_replace_all(., "__", "_")) %>% # remove double "_"
    column_to_rownames()
  )
    
  ALDEx2plot_pairwise <- CalcALDEx(
    physeq_obj = Ps_obj_filt_median_tea_pairwise_s,
    vars2test = "Season",
    rare_phyla = Rare_phyla,
    sig_level = significance,
    LFC = 0
    )
  ALDEx_comparisons$Results[[j]] <- ALDEx2plot_pairwise # store results
  ALDEx_comparisons$Results[[j]]$Var1 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[1]
  ALDEx_comparisons$Results[[j]]$Var2 <- str_split(ALDEx_comparisons$Comparisons[j], " vs ", simplify = TRUE)[2]
  ALDEX_summary <- tibble(Label = c(paste0("⬆", sum( ALDEx_comparisons$Results[[j]]$effect > 0 &  ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " ⬇", sum( ALDEx_comparisons$Results[[j]]$effect < 0 &  ALDEx_comparisons$Results[[j]]$Significance == "Pass"), " (", nrow( ALDEx_comparisons$Results[[j]]), ")")))
  
  ALDEx2plot_pairwise %>%
    filter(Significance == "Pass") %>%
    dplyr::select(OTU, baseMean, effect, Phylum, Class, Order, Family, Genus) %>%
    arrange(desc(abs(effect))) ->
    ALDEx2plot_pairwise_results
    
  # print(ALDEx2plot_pairwise_results %>% 
  #   kable(., digits = c(2), caption = "Significantly different taxa:") %>%
  #   kable_styling(
  #     bootstrap_options = c("striped", "hover", "condensed", "responsive"),
  #     full_width = F
  #   ))
  
  write.csv(ALDEx2plot_pairwise_results, file = paste0("Aldex", "_", paste0(comparison_string, collapse = "_"), ".csv"))
  
  # Plot ALDEX plot
  p1 <- GGPlotALDExTax(ALDEx2plot_pairwise, OTU_labels = FALSE, sig_level = significance) +
    # ggtitle(ALDEx_comparisons$Comparisons[j]) +
    geom_text(
    data    = ALDEX_summary,
    mapping = aes(x = Inf, y = Inf, label = Label),
    hjust   = 1.1,
    vjust   = 1.6
  ) 
  p1 <- p1 + labs(title = ALDEx_comparisons$Comparisons[j])
  print(p1)
  
  # Plot OTU plots
  p2 <- GGPlotTopOTUs(
    physeq_obj = Ps_obj_filt_median_tea_pairwise_s,
    vars2test = "Season",
    ALDEx_obj = ALDEx2plot_pairwise,
    rank_by = "effect",
    Ntop = 12
  )
  print(p2)
}

phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 2700 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 2700 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 2700 tips and 2698 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 45 samples ] sample_data() Sample Data: [ 45 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3875 taxa and 45 samples ] sample_data() Sample Data: [ 45 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3875 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3875 tips and 3873 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 48 samples ] sample_data() Sample Data: [ 48 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 3270 taxa and 48 samples ] sample_data() Sample Data: [ 48 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 3270 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 3270 tips and 3268 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 44 samples ] sample_data() Sample Data: [ 44 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 2511 taxa and 44 samples ] sample_data() Sample Data: [ 44 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 2511 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 2511 tips and 2509 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 1818 taxa and 47 samples ] sample_data() Sample Data: [ 47 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 1818 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 1818 tips and 1816 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 4217 taxa and 45 samples ] sample_data() Sample Data: [ 45 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 4217 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 4217 tips and 4215 internal nodes ] phyloseq-class experiment-level object otu_table() OTU Table: [ 2745 taxa and 45 samples ] sample_data() Sample Data: [ 45 samples by 25 sample variables ] tax_table() Taxonomy Table: [ 2745 taxa by 6 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 2745 tips and 2743 internal nodes ]

filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU[common_taxa <- filter(ALDEx_comparisons$Results[[1]], Significance == "Pass")$OTU %in% filter(ALDEx_comparisons$Results[[2]], Significance == "Pass")$OTU]

[1] "Seq_511" "Seq_35" "Seq_99" "Seq_119" "Seq_540" "Seq_87" "Seq_477" [8] "Seq_1578" "Seq_169" "Seq_371" "Seq_216" "Seq_983" "Seq_170" "Seq_1351" [15] "Seq_526" "Seq_536" "Seq_644" "Seq_2982" "Seq_1503" "Seq_1904" "Seq_2209" [22] "Seq_222" "Seq_193" "Seq_2377" "Seq_2325" "Seq_552" "Seq_1633" "Seq_321" [29] "Seq_1632" "Seq_550" "Seq_796" "Seq_297" "Seq_1548" "Seq_1581" "Seq_1866" [36] "Seq_1658" "Seq_265" "Seq_965" "Seq_1569" "Seq_151" "Seq_2024" "Seq_2469" [43] "Seq_276" "Seq_1539" "Seq_2269" "Seq_667" "Seq_381" "Seq_340" "Seq_584" [50] "Seq_211" "Seq_365" "Seq_661" "Seq_734" "Seq_615" "Seq_416" "Seq_81"
[57] "Seq_153" "Seq_1669" "Seq_1264" "Seq_1432" "Seq_613" "Seq_1353" "Seq_217" [64] "Seq_1911" "Seq_66" "Seq_791" "Seq_1895" "Seq_780" "Seq_201" "Seq_527" [71] "Seq_267" "Seq_47" "Seq_94" "Seq_9" "Seq_275" "Seq_424" "Seq_741" [78] "Seq_1731" "Seq_606" "Seq_113" "Seq_55" "Seq_16" "Seq_134" "Seq_2503" [85] "Seq_1173" "Seq_29" "Seq_637" "Seq_18" "Seq_27" "Seq_232" "Seq_101" [92] "Seq_608" "Seq_164" "Seq_1326" "Seq_1025" "Seq_1704" "Seq_481" "Seq_149" [99] "Seq_619" "Seq_1118" "Seq_1525" "Seq_486" "Seq_96" "Seq_1024" "Seq_50"
[106] "Seq_286" "Seq_60" "Seq_618" "Seq_364" "Seq_647" "Seq_531" "Seq_790" [113] "Seq_1267" "Seq_2041" "Seq_2228" "Seq_1961" "Seq_1197" "Seq_2482" "Seq_187" [120] "Seq_336" "Seq_979" "Seq_70" "Seq_285" "Seq_175" "Seq_991" "Seq_2045" [127] "Seq_713" "Seq_1108" "Seq_927" "Seq_625" "Seq_1255" "Seq_263" "Seq_582" [134] "Seq_249" "Seq_108" "Seq_832" "Seq_366" "Seq_329" "Seq_215" "Seq_376" [141] "Seq_464" "Seq_567" "Seq_32" "Seq_442" "Seq_1320" "Seq_1471" "Seq_764" [148] "Seq_851" "Seq_1732" "Seq_385" "Seq_257" "Seq_161" "Seq_341" "Seq_107" [155] "Seq_139" "Seq_1071" "Seq_17" "Seq_15" "Seq_54" "Seq_30" "Seq_82"
[162] "Seq_431" "Seq_74" "Seq_1354" "Seq_1579" "Seq_157" "Seq_330" "Seq_71"
[169] "Seq_98" "Seq_459" "Seq_7" "Seq_1613" "Seq_441" "Seq_537" "Seq_1750" [176] "Seq_704" "Seq_752" "Seq_465" "Seq_444" "Seq_312" "Seq_466" "Seq_669" [183] "Seq_224" "Seq_214" "Seq_358" "Seq_600" "Seq_627" "Seq_427" "Seq_180" [190] "Seq_351" "Seq_200" "Seq_210" "Seq_142" "Seq_1903" "Seq_273" "Seq_1226" [197] "Seq_1387" "Seq_127" "Seq_255" "Seq_57" "Seq_191" "Seq_1286" "Seq_199" [204] "Seq_1477" "Seq_13" "Seq_1681" "Seq_42" "Seq_724" "Seq_229" "Seq_33"
[211] "Seq_598" "Seq_414" "Seq_296" "Seq_2011" "Seq_2170" "Seq_2004" "Seq_357" [218] "Seq_189" "Seq_91" "Seq_1113" "Seq_413" "Seq_344" "Seq_575" "Seq_1980" [225] "Seq_360" "Seq_59" "Seq_828" "Seq_2014" "Seq_242" "Seq_2856" "Seq_793" [232] "Seq_31" "Seq_41" "Seq_168" "Seq_230" "Seq_293" "Seq_1053" "Seq_877" [239] "Seq_36" "Seq_39" "Seq_62" "Seq_93" "Seq_205" "Seq_394" "Seq_103" [246] "Seq_73" "Seq_162" "Seq_1997" "Seq_684" "Seq_546" "Seq_530" "Seq_1020" [253] "Seq_1175" "Seq_1257" "Seq_693" "Seq_484" "Seq_642" "Seq_1473" "Seq_188" [260] "Seq_1612" "Seq_288" "Seq_208" "Seq_646" "Seq_2110" "Seq_271" "Seq_184" [267] "Seq_1317" "Seq_318" "Seq_136" "Seq_48" "Seq_146" "Seq_95" "Seq_523" [274] "Seq_158" "Seq_1134" "Seq_811" "Seq_90" "Seq_156" "Seq_308" "Seq_290" [281] "Seq_788" "Seq_2983" "Seq_110" "Seq_1915" "Seq_58" "Seq_1360" "Seq_1050" [288] "Seq_1090" "Seq_1135" "Seq_471" "Seq_369" "Seq_1033" "Seq_11" "Seq_458" [295] "Seq_595" "Seq_2638" "Seq_114"

sessioninfo::session_info() %>%
  details::details(
    summary = 'Current session info',
    open    = TRUE
 )

Current session info


─ Session info ─────────────────────────────────────────────────────────────────────────
 setting  value                       
 version  R version 4.0.3 (2020-10-10)
 os       Ubuntu 18.04.5 LTS          
 system   x86_64, linux-gnu           
 ui       X11                         
 language en_GB                       
 collate  en_DK.UTF-8                 
 ctype    en_DK.UTF-8                 
 tz       Europe/Prague               
 date     2021-05-21                  

─ Packages ─────────────────────────────────────────────────────────────────────────────
 package              * version  date       lib source                                  
 ade4                 * 1.7-16   2020-10-28 [1] CRAN (R 4.0.2)                          
 affy                   1.68.0   2020-10-27 [1] Bioconductor                            
 affyio                 1.60.0   2020-10-27 [1] Bioconductor                            
 ALDEx2               * 1.22.0   2020-10-27 [1] Bioconductor                            
 ape                    5.5      2021-04-25 [1] CRAN (R 4.0.3)                          
 assertthat             0.2.1    2019-03-21 [1] CRAN (R 4.0.2)                          
 backports              1.2.1    2020-12-09 [1] CRAN (R 4.0.2)                          
 Biobase              * 2.50.0   2020-10-27 [1] Bioconductor                            
 BiocGenerics         * 0.36.1   2021-04-16 [1] Bioconductor                            
 BiocManager            1.30.15  2021-05-11 [1] CRAN (R 4.0.3)                          
 BiocParallel           1.24.1   2020-11-06 [1] Bioconductor                            
 biomformat             1.18.0   2020-10-27 [1] Bioconductor                            
 Biostrings             2.58.0   2020-10-27 [1] Bioconductor                            
 bitops                 1.0-7    2021-04-24 [1] CRAN (R 4.0.3)                          
 broom                * 0.7.6    2021-04-05 [1] CRAN (R 4.0.3)                          
 bslib                  0.2.5.1  2021-05-18 [1] CRAN (R 4.0.3)                          
 cellranger             1.1.0    2016-07-27 [1] CRAN (R 4.0.2)                          
 cli                    2.5.0    2021-04-26 [1] CRAN (R 4.0.3)                          
 clipr                  0.7.1    2020-10-08 [1] CRAN (R 4.0.2)                          
 cluster                2.1.2    2021-04-17 [1] CRAN (R 4.0.3)                          
 codetools              0.2-18   2020-11-04 [1] CRAN (R 4.0.2)                          
 colorspace             2.0-1    2021-05-04 [1] CRAN (R 4.0.3)                          
 cowplot              * 1.1.1    2020-12-30 [1] CRAN (R 4.0.2)                          
 crayon                 1.4.1    2021-02-08 [1] CRAN (R 4.0.3)                          
 data.table             1.14.0   2021-02-21 [1] CRAN (R 4.0.3)                          
 DBI                    1.1.1    2021-01-15 [1] CRAN (R 4.0.3)                          
 dbplyr                 2.1.1    2021-04-06 [1] CRAN (R 4.0.3)                          
 DelayedArray           0.16.3   2021-03-24 [1] Bioconductor                            
 desc                   1.3.0    2021-03-05 [1] CRAN (R 4.0.3)                          
 details                0.2.1    2020-01-12 [1] CRAN (R 4.0.2)                          
 digest                 0.6.27   2020-10-24 [1] CRAN (R 4.0.2)                          
 doParallel           * 1.0.16   2020-10-16 [1] CRAN (R 4.0.2)                          
 dplyr                * 1.0.6    2021-05-05 [1] CRAN (R 4.0.3)                          
 ellipsis               0.3.2    2021-04-29 [1] CRAN (R 4.0.3)                          
 evaluate               0.14     2019-05-28 [1] CRAN (R 4.0.2)                          
 extrafont            * 0.17     2014-12-08 [1] CRAN (R 4.0.2)                          
 extrafontdb            1.0      2012-06-11 [1] CRAN (R 4.0.2)                          
 fansi                  0.4.2    2021-01-15 [1] CRAN (R 4.0.3)                          
 farver                 2.1.0    2021-02-28 [1] CRAN (R 4.0.3)                          
 forcats              * 0.5.1    2021-01-27 [1] CRAN (R 4.0.3)                          
 foreach              * 1.5.1    2020-10-15 [1] CRAN (R 4.0.2)                          
 fs                     1.5.0    2020-07-31 [1] CRAN (R 4.0.2)                          
 generics               0.1.0    2020-10-31 [1] CRAN (R 4.0.2)                          
 GenomeInfoDb           1.26.7   2021-04-08 [1] Bioconductor                            
 GenomeInfoDbData       1.2.4    2021-05-05 [1] Bioconductor                            
 GenomicRanges          1.42.0   2020-10-27 [1] Bioconductor                            
 ggplot2              * 3.3.3    2020-12-30 [1] CRAN (R 4.0.2)                          
 ggpomological        * 0.1.2    2020-08-13 [1] Github (gadenbuie/ggpomological@69f3815)
 ggrepel              * 0.9.1    2021-01-15 [1] CRAN (R 4.0.3)                          
 ggsci                * 2.9      2018-05-14 [1] CRAN (R 4.0.2)                          
 ggthemes             * 4.2.4    2021-01-20 [1] CRAN (R 4.0.3)                          
 glue                   1.4.2    2020-08-27 [1] CRAN (R 4.0.2)                          
 gridExtra            * 2.3      2017-09-09 [1] CRAN (R 4.0.2)                          
 gtable                 0.3.0    2019-03-25 [1] CRAN (R 4.0.2)                          
 haven                  2.4.1    2021-04-23 [1] CRAN (R 4.0.3)                          
 hexbin                 1.28.2   2021-01-08 [1] CRAN (R 4.0.2)                          
 highr                  0.9      2021-04-16 [1] CRAN (R 4.0.3)                          
 hms                    1.1.0    2021-05-17 [1] CRAN (R 4.0.3)                          
 htmltools              0.5.1.1  2021-01-22 [1] CRAN (R 4.0.3)                          
 httr                   1.4.2    2020-07-20 [1] CRAN (R 4.0.2)                          
 igraph                 1.2.6    2020-10-06 [1] CRAN (R 4.0.2)                          
 IRanges                2.24.1   2020-12-12 [1] Bioconductor                            
 iterators            * 1.0.13   2020-10-15 [1] CRAN (R 4.0.2)                          
 jquerylib              0.1.4    2021-04-26 [1] CRAN (R 4.0.3)                          
 jsonlite               1.7.2    2020-12-09 [1] CRAN (R 4.0.2)                          
 kableExtra           * 1.3.4    2021-02-20 [1] CRAN (R 4.0.3)                          
 knitr                * 1.33     2021-04-24 [1] CRAN (R 4.0.3)                          
 labeling               0.4.2    2020-10-20 [1] CRAN (R 4.0.2)                          
 lattice              * 0.20-44  2021-05-02 [1] CRAN (R 4.0.3)                          
 lifecycle              1.0.0    2021-02-15 [1] CRAN (R 4.0.3)                          
 limma                  3.46.0   2020-10-27 [1] Bioconductor                            
 lubridate              1.7.10   2021-02-26 [1] CRAN (R 4.0.3)                          
 magrittr             * 2.0.1    2020-11-17 [1] CRAN (R 4.0.2)                          
 MASS                 * 7.3-54   2021-05-03 [1] CRAN (R 4.0.3)                          
 Matrix                 1.3-3    2021-05-04 [1] CRAN (R 4.0.3)                          
 MatrixGenerics         1.2.1    2021-01-30 [1] Bioconductor                            
 matrixStats            0.58.0   2021-01-29 [1] CRAN (R 4.0.3)                          
 mgcv                   1.8-35   2021-04-18 [1] CRAN (R 4.0.3)                          
 modelr                 0.1.8    2020-05-19 [1] CRAN (R 4.0.2)                          
 multtest               2.46.0   2020-10-27 [1] Bioconductor                            
 munsell                0.5.0    2018-06-12 [1] CRAN (R 4.0.2)                          
 NADA                 * 1.6-1.1  2020-03-22 [1] CRAN (R 4.0.2)                          
 nlme                   3.1-152  2021-02-04 [1] CRAN (R 4.0.3)                          
 permute              * 0.9-5    2019-03-12 [1] CRAN (R 4.0.2)                          
 phyloseq             * 1.34.0   2020-10-27 [1] Bioconductor                            
 pillar                 1.6.1    2021-05-16 [1] CRAN (R 4.0.3)                          
 pkgconfig              2.0.3    2019-09-22 [1] CRAN (R 4.0.2)                          
 plyr                   1.8.6    2020-03-03 [1] CRAN (R 4.0.2)                          
 png                    0.1-7    2013-12-03 [1] CRAN (R 4.0.2)                          
 preprocessCore         1.52.1   2021-01-08 [1] Bioconductor                            
 prettyunits            1.1.1    2020-01-24 [1] CRAN (R 4.0.2)                          
 progress               1.2.2    2019-05-16 [1] CRAN (R 4.0.2)                          
 ps                     1.6.0    2021-02-28 [1] CRAN (R 4.0.3)                          
 purrr                * 0.3.4    2020-04-17 [1] CRAN (R 4.0.2)                          
 R6                     2.5.0    2020-10-28 [1] CRAN (R 4.0.2)                          
 Rcpp                   1.0.6    2021-01-15 [1] CRAN (R 4.0.3)                          
 RCurl                  1.98-1.3 2021-03-16 [1] CRAN (R 4.0.3)                          
 readr                * 1.4.0    2020-10-05 [1] CRAN (R 4.0.2)                          
 readxl                 1.3.1    2019-03-13 [1] CRAN (R 4.0.2)                          
 reprex                 2.0.0    2021-04-02 [1] CRAN (R 4.0.3)                          
 reshape2               1.4.4    2020-04-09 [1] CRAN (R 4.0.2)                          
 rhdf5                  2.34.0   2020-10-27 [1] Bioconductor                            
 rhdf5filters           1.2.1    2021-05-03 [1] Bioconductor                            
 Rhdf5lib               1.12.1   2021-01-26 [1] Bioconductor                            
 rlang                  0.4.11   2021-04-30 [1] CRAN (R 4.0.3)                          
 rmarkdown            * 2.8      2021-05-07 [1] CRAN (R 4.0.3)                          
 rprojroot              2.0.2    2020-11-15 [1] CRAN (R 4.0.2)                          
 rstudioapi             0.13     2020-11-12 [1] CRAN (R 4.0.2)                          
 Rttf2pt1               1.3.8    2020-01-10 [1] CRAN (R 4.0.2)                          
 rvest                  1.0.0    2021-03-09 [1] CRAN (R 4.0.3)                          
 S4Vectors              0.28.1   2020-12-09 [1] Bioconductor                            
 sass                   0.4.0    2021-05-12 [1] CRAN (R 4.0.3)                          
 scales               * 1.1.1    2020-05-11 [1] CRAN (R 4.0.2)                          
 sessioninfo            1.1.1    2018-11-05 [1] CRAN (R 4.0.2)                          
 stringi                1.6.2    2021-05-17 [1] CRAN (R 4.0.3)                          
 stringr              * 1.4.0    2019-02-10 [1] CRAN (R 4.0.2)                          
 SummarizedExperiment   1.20.0   2020-10-27 [1] Bioconductor                            
 survival             * 3.2-11   2021-04-26 [1] CRAN (R 4.0.3)                          
 svglite              * 2.0.0    2021-02-20 [1] CRAN (R 4.0.3)                          
 systemfonts            1.0.2    2021-05-11 [1] CRAN (R 4.0.3)                          
 tibble               * 3.1.2    2021-05-16 [1] CRAN (R 4.0.3)                          
 tidyr                * 1.1.3    2021-03-03 [1] CRAN (R 4.0.3)                          
 tidyselect             1.1.1    2021-04-30 [1] CRAN (R 4.0.3)                          
 tidyverse            * 1.3.1    2021-04-15 [1] CRAN (R 4.0.3)                          
 truncnorm            * 1.0-8    2018-02-27 [1] CRAN (R 4.0.2)                          
 utf8                   1.2.1    2021-03-12 [1] CRAN (R 4.0.3)                          
 vctrs                  0.3.8    2021-04-29 [1] CRAN (R 4.0.3)                          
 vegan                * 2.5-7    2020-11-28 [1] CRAN (R 4.0.3)                          
 viridisLite            0.4.0    2021-04-13 [1] CRAN (R 4.0.3)                          
 vsn                  * 3.58.0   2020-10-27 [1] Bioconductor                            
 webshot                0.5.2    2019-11-22 [1] CRAN (R 4.0.2)                          
 withr                  2.4.2    2021-04-18 [1] CRAN (R 4.0.3)                          
 xfun                   0.23     2021-05-15 [1] CRAN (R 4.0.3)                          
 xml2                   1.3.2    2020-04-23 [1] CRAN (R 4.0.2)                          
 XVector                0.30.0   2020-10-27 [1] Bioconductor                            
 yaml                   2.2.1    2020-02-01 [1] CRAN (R 4.0.2)                          
 zCompositions        * 1.3.4    2020-03-04 [1] CRAN (R 4.0.2)                          
 zlibbioc               1.36.0   2020-10-27 [1] Bioconductor                            

[1] /home/angel/R/library
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library


References

Fernandes AD, Macklaim JM, Linn TG et al. ANOVA-Like Differential Expression (ALDEx) Analysis for Mixed Population RNA-Seq. PLOS ONE 2013;8:e67019.